### Calculate counterfactuals related to military coalitions from the base model

date()
library("tidyverse")
library("assertr")
library("foreach")
library("geosphere")
library("RColorBrewer")
source("backend_main.r")
sessionInfo()


load("results/fit_base.rda")
load("kr_dyad_year.rda")
load("kr_mid_distance.rda")
load("kr_state_year.rda")


## -----------------------------------------------------------------------------
## General functions
## -----------------------------------------------------------------------------

## Create "synthetic" dispute participant data from the imputed state-year data
synth_participant <- foreach (imp = imp_state_year) %do% {
  mutate(imp, sidea = NA, id = year)
}

avg_from_list <- function(x, which) {
  x <- map(x, ~ pluck(., which))
  n <- length(x)
  x <- reduce(x, ~ .x + .y)
  x / n
}

coalition_cf <- function(dispute_dyads,
                         fit_list,
                         mid_id = NULL,
                         latlong = NULL,
                         comment = NULL,
                         synth_participant_ = synth_participant) {
  ## Extract each state's side and first year in dispute
  ##
  ## n.b., states only switch sides in WW2 in our data
  participants <- dispute_dyads %>%
    select(ccode = ccode1, year, sidea = sidea1) %>%
    rbind(dispute_dyads %>% select(ccode = ccode2, year, sidea = sidea2)) %>%
    add_column(id = 1, .before = 1) %>%
    distinct() %>%
    assert(not_na, everything())

  ## Take a list of dyad dispute-years (typically the output of
  ## make_dispute_dyads) and generate dispute data mirroring the data used to
  ## fit the model
  newdata_dispute <- foreach (imp_d = imp_dyad_year, imp_s = imp_state_year) %do% {
    ## Calculate dispute-level variables pertaining to each side
    vars_by_side <- participants %>%
      left_join(imp_s, by = c("ccode", "year")) %>%
      summarise(polity_a = mean(polity2[sidea == 1]),
                polity_b = mean(polity2[sidea == 0]),
                majpow_a = max(majpow[sidea == 1]),
                majpow_b = max(majpow[sidea == 0]),
                n_states_a = sum(sidea == 1),
                n_states_b = sum(sidea == 0))

    ## Calculate other dispute-level variables and merge in the side-specific
    ## ones calculated above
    dispute_dyads %>%
      select(ccode1, ccode2, year) %>%
      assert(not_na, everything()) %>%
      verify(ccode1 < ccode2) %>%
      select(ccode_a = ccode1,
             ccode_b = ccode2,
             year) %>%
      left_join(imp_d, by = c("ccode_a", "ccode_b", "year")) %>%
      assert(not_na, everything()) %>%
      verify(!duplicated(paste(ccode_a, ccode_b))) %>%
      summarise(py_alt = mean(py_alt),
                s_cinc = mean(s_cinc),
                contig = mean(contig)) %>%
      bind_cols(vars_by_side) %>%
      add_column(id = 1, .before = 1) %>%
      add_column(war = 1) %>%
      add_column(win_a = 0) %>%
      add_column(win_b = 0) %>%
      add_column(win_b_alt = 0)
  }

  ## Calculate distance to the dispute
  if (!is.null(mid_id)) {
    dat_dist <- data_mid_distance %>%
      filter(id == !!mid_id) %>%
      select(ccode, distance)
  } else {
    dat_dist <- participants %>%
      left_join(imp_state_year[[1]] %>% select(ccode, year, latitude, longitude),
                by = c("ccode", "year")) %>%
      mutate(distance = distGeo(cbind(longitude, latitude),
                                cbind(!!latlong[2], !!latlong[1])),
             distance = distance / 1609.344) %>%
      select(ccode, distance)
  }

  ## Merge in participant-level variables as well as distance to the dispute
  ## (as calculated in our original data)
  newdata_participant <- foreach (imp = synth_participant_) %do% {
    participants %>%
      left_join(imp %>% select(-id, -stateabb, -statenme, -sidea),
                by = c("ccode", "year")) %>%
      left_join(dat_dist, by = c("ccode")) %>%
      assert(not_na, -ends_with("_lag_pure")) %>%
      verify(!duplicated(paste(ccode, year)))
  }

  ans <- foreach (fit = fit_list, nd_disp = newdata_dispute, nd_part = newdata_participant) %do% {
    ## Extract model matrices
    setup <- structwar_setup(f_dispute = fit[["f_dispute"]],
                             f_participant = fit[["f_participant"]],
                             data_dispute = nd_disp,
                             data_participant = nd_part,
                             xlev_dispute = fit[["xlev_dispute"]],
                             xlev_participant = fit[["xlev_participant"]],
                             n_halton = fit[["n_halton"]])
    X <- setup[["X"]]
    Z <- setup[["Z"]]
    L_a <- setup[["L_a"]]
    L_b <- setup[["L_b"]]
    S_a <- setup[["S_a"]]
    S_b <- setup[["S_b"]]

    ## Calculate model parameters
    cf <- coef(fit[["fit"]])
    parms <- extract_params(est = cf,
                            setup = setup,
                            for_counterfactuals = TRUE)
    beta <- cf[seq_len(ncol(X))]
    gamma <- cf[ncol(X) + seq_len(ncol(Z))]

    ## Calculate Nash equilibrium of the contest subgame
    ratio <- parms$ratio
    ratio_sort <- sort(ratio)
    criterion <- cumsum(ratio_sort) - (seq_along(ratio_sort) - 1) * ratio_sort
    J <- sum(criterion > 0)
    stopifnot(J >= 2)
    max_ratio <- ratio_sort[J]
    active <- !(ratio > max_ratio)
    ratio_denom <- sum(ratio[active])
    q <- (J - 1) * ratio / ratio_denom
    cc <- parms$cost
    e <- if_else(active, q * (1 - q) / cc, 0)
    p <- if_else(active, 1 - q, 0)
    pi <- p^2

    ## Calculate bargaining equilibrium
    sidea <- participants$sidea
    pi_a <- sum(pi[sidea == 1])
    pi_b <- sum(pi[sidea == 0])
    pr_win_a <- sum(p[sidea == 1])
    pwax <- prob_war_and_x(p_theta_a = setup$p_theta_a,
                           loc_a = parms$loc_a,
                           loc_b = parms$loc_b,
                           scl_a = parms$scl_a,
                           scl_b = parms$scl_b,
                           pi_a = pi_a,
                           pi_b = pi_b)
    x_star <- pwax$x_star
    pr_war <- pwax$pr_war

    ans_participant <- participants %>%
      add_column(cost = cc, ratio = ratio, effort = e, wt_effort = cc * e / ratio, p, pi) %>%
      select(-id)

    ans_dispute <-
      tibble(pr_win_a, x_star, pr_war, loc_a = parms$loc_a, loc_b = parms$loc_b,
             scl_a = parms$scl_a, scl_b = parms$scl_b)

    list(participant = ans_participant,
         dispute = ans_dispute,
         X = X,
         Z = Z,
         beta = beta,
         gamma = gamma)
  }

  ans_participant <- lapply(ans, "[[", "participant")
  ans_participant <- do.call(rbind, ans_participant) %>%
    group_by(ccode) %>%
    summarise_all(~ mean(.))

  ans_dispute <- lapply(ans, "[[", "dispute")
  ans_dispute <- do.call(rbind, ans_dispute) %>%
    summarise_all(~ mean(.))

  if (!is.null(comment)) {
    ans_participant <- add_column(ans_participant, comment = !! comment)
    ans_dispute <- add_column(ans_dispute, comment = !! comment)
  }

  list(participant = ans_participant,
       dispute = ans_dispute,
       X = avg_from_list(ans, "X"),
       Z = avg_from_list(ans, "Z"),
       beta = avg_from_list(ans, "beta"),
       gamma = avg_from_list(ans, "gamma"))
}

make_dispute_dyads <- function(a, b, year) {
  crossing(a, b, year) %>%
    mutate(ccode1 = pmin(a, b),
           ccode2 = pmax(a, b),
           sidea1 = as.numeric(ccode1 %in% !! a),
           sidea2 = as.numeric(ccode2 %in% !! a)) %>%
    verify(sidea1 + sidea2 == 1)
}

collapse_results <- function(...) {
  dots <- list(...)
  ans_part <- do.call(rbind, lapply(dots, "[[", "participant"))
  ans_disp <- do.call(rbind, lapply(dots, "[[", "dispute"))
  list(participant = ans_part, dispute = ans_disp)
}

make_table <- function(dat, win_name = "") {
  dat <- dat %>%
    mutate(pr_win_a = sprintf("%.2f", pr_win_a),
           pr_war = sprintf("%.2f", pr_war),
           entry = str_glue("{comment} & {pr_win_a} \\\\"))
  c("\\begin{tabular}{lcc}",
    "\\toprule",
    str_glue("Counterfactual & Pr({win_name}) \\\\"),
    "\\midrule",
    dat$entry,
    "\\bottomrule",
    "\\end{tabular}")
}

## Convenience function for comparing two countries' results from a
## counterfactual object
compare <- function(cf, ccode1, ccode2) {
  idx_1 <- which(cf[["participant"]]$ccode == ccode1)
  idx_2 <- which(cf[["participant"]]$ccode == ccode2)
  x_1 <- drop(cbind(cf$X[idx_1, , drop = FALSE], cf$Z[idx_1, , drop = FALSE]))
  x_2 <- drop(cbind(cf$X[idx_2, , drop = FALSE], cf$Z[idx_2, , drop = FALSE]))
  beta <- c(cf$beta, -1 * cf$gamma)
  tibble(term = !! names(beta),
         beta = !! beta,
         !! paste0("x_", ccode1) := !! x_1,
         !! paste0("x_", ccode2) := !! x_2,
         x_diff = !! (x_1 - x_2),
         !! paste0("x_", ccode1, "_wt") := !! x_1 * beta,
         !! paste0("x_", ccode2, "_wt") := !! x_2 * beta,
         x_diff_wt = !! ((x_1 - x_2) * beta))
}

country_names <- tribble(
  ~ccode, ~country,
  2, "USA",
  200, "UK",
  210, "Netherlands",
  211, "Belgium",
  220, "France",
  230, "Spain",
  255, "Germany",
  290, "Poland",
  300, "Austria",
  325, "Italy",
  345, "Serbia",
  350, "Greece",
  355, "Bulgaria",
  360, "Romania",
  365, "Russia",
  390, "Denmark",
  640, "Turkey",
  645, "Iraq",
  710, "China",
  713, "Taiwan",
  740, "Japan",
  900, "Australia"
)


## -----------------------------------------------------------------------------
## World War I (257)
## -----------------------------------------------------------------------------

dat_ww1_full <- make_dispute_dyads(a = c(2, 200, 211, 220, 345, 365),
                                   b = c(255, 300, 640),
                                   year = 1914)
cf_ww1_full <- coalition_cf(dispute_dyads = dat_ww1_full,
                            fit_list = fit_base,
                            mid_id = 257,
                            comment = "With USA")

dat_ww1_no_usa <- dat_ww1_full %>%
  filter(ccode1 != 2, ccode2 != 2)
cf_ww1_no_usa <- coalition_cf(dispute_dyads = dat_ww1_no_usa,
                              fit_list = fit_base,
                              mid_id = 257,
                              comment = "Baseline")

dat_ww1_no_rus <- dat_ww1_full %>%
  filter(!(ccode1 %in% c(2, 365)), !(ccode2 %in% c(2, 365)))
cf_ww1_no_rus <- coalition_cf(dispute_dyads = dat_ww1_no_rus,
                              fit_list = fit_base,
                              mid_id = 257,
                              comment = "Without Russia")

dat_ww1_no_ukg <- dat_ww1_full %>%
  filter(!(ccode1 %in% c(2, 200)), !(ccode2 %in% c(2, 200)))
cf_ww1_no_ukg <- coalition_cf(dispute_dyads = dat_ww1_no_ukg,
                              fit_list = fit_base,
                              mid_id = 257,
                              comment = "Without UK")

dat_ww1 <- collapse_results(cf_ww1_no_usa, cf_ww1_no_rus, cf_ww1_no_ukg, cf_ww1_full)
dat_ww1_part <- dat_ww1[["participant"]]
dat_ww1_disp <- dat_ww1[["dispute"]]

pdat_ww1_part <-
  complete(dat_ww1_part, ccode, comment, fill = list(effort = 0, wt_effort = 0)) %>%
  left_join(country_names, by = "ccode") %>%
  mutate(alliance = if_else(ccode %in% c(255, 300, 640), "Central Powers", "Allied Powers"),
         Players = factor(comment, levels = c("Baseline", "Without UK", "Without Russia", "With USA")))

plot_ww1 <- ggplot(pdat_ww1_part, aes(x = country, y = wt_effort)) +
  geom_col(aes(fill = Players),
           position = position_dodge2(preserve = "single")) +
  facet_grid(~ alliance, scales = "free_x", space = "free") +
  scale_x_discrete("") +
  scale_y_continuous("Equilibrium contribution") +
  scale_fill_manual(values = brewer.pal(5, "Reds")[-1]) +
  ggtitle("Counterfactuals for World War I, 1914")

if (!dir.exists("figures"))
  dir.create("figures")
pdf(file = "figures/figure_2.pdf",
    width = 7,
    height = 3.25)
print(plot_ww1)
dev.off()

if (!dir.exists("tables"))
  dir.create("tables")
tab_ww1 <- make_table(dat_ww1_disp, win_name = "Allied Powers Win")
writeLines(tab_ww1, con = "tables/table_5a.tex")


## -----------------------------------------------------------------------------
## World War I, if Schlieffen Plan had succeeded in 1914
## -----------------------------------------------------------------------------

dat_sch_base <- make_dispute_dyads(a = c(200, 211, 220, 345, 365),
                                   b = c(255, 300, 640),
                                   year = 1915)
cf_sch_base <- coalition_cf(dispute_dyads = dat_sch_base,
                            fit_list = fit_base,
                            mid_id = 257,
                            comment = "Baseline")

dat_sch_no_bel <- dat_sch_base %>%
  filter(ccode1 != 211, ccode2 != 211)
cf_sch_no_bel <- coalition_cf(dispute_dyads = dat_sch_no_bel,
                              fit_list = fit_base,
                              mid_id = 257,
                              comment = "Without Belgium")

dat_sch_no_fr <- dat_sch_base %>%
  filter(ccode1 != 220, ccode2 != 220)
cf_sch_no_fr <- coalition_cf(dispute_dyads = dat_sch_no_fr,
                             fit_list = fit_base,
                             mid_id = 257,
                             comment = "Without France")

dat_sch_no_both <- inner_join(dat_sch_no_bel, dat_sch_no_fr)
cf_sch_no_both <- coalition_cf(dispute_dyads = dat_sch_no_both,
                               fit_list = fit_base,
                               mid_id = 257,
                               comment = "Without both")

## Calculate USA utility on entry in April 1917
##
## Changes from 1914 composition by now:
##   Italy joined Entente April 1915
##   Bulgaria joined Triple Alliance October 1915
##   Serbia defeated November 1915
##   Romania joined Entente August 1916
dat_ww1_1917 <- make_dispute_dyads(a = c(2, 200, 211, 220, 325, 360, 365),
                                   b = c(255, 300, 355, 640),
                                   year = 1917)
cf_ww1_1917 <- coalition_cf(dispute_dyads = dat_ww1_1917,
                                      fit_list = fit_base,
                                      mid_id = 257)

## Calculate USA utility from hypothetical entry if France and Belgium had been knocked out by the
## Schlieffen plan
dat_sch_usa <- make_dispute_dyads(a = c(2, 200, 345, 365),
                                  b = c(255, 300, 640),
                                  year = 1915)
cf_sch_usa <- coalition_cf(dispute_dyads = dat_sch_usa,
                                     fit_list = fit_base,
                                     mid_id = 257,
                                     comment = "Without both, with USA")

cat("\nUSA utility from entering WW1 in 1917:\n")
print(cf_ww1_1917$participant %>% filter(ccode == 2) %>% pull(pi))
cat("\nUSA utility from entering WW1 in 1915 if Schlieffen Plan worked:\n")
print(cf_sch_usa$participant %>% filter(ccode == 2) %>% pull(pi))
cat("\nEntente probability of victory if US entered:\n")
print(cf_sch_usa$dispute)

dat_sch <- collapse_results(cf_sch_base, cf_sch_no_bel, cf_sch_no_fr, cf_sch_no_both, cf_sch_usa)
dat_sch_part <- dat_sch[["participant"]]
dat_sch_disp <- dat_sch[["dispute"]]

tab_sch <- make_table(dat_sch_disp, win_name = "Allied Powers Win")
writeLines(tab_sch, con = "tables/table_5b.tex")


## -----------------------------------------------------------------------------
## Iraq war crisis (bargaining counterfactuals)
## -----------------------------------------------------------------------------

willing_coalition <- c(2, 200, 290, 390, 900)

## Just USA and Iraq
dat_iraq_bilat <- make_dispute_dyads(a = 2, b = 645, year = 2003)
latlong_iraq <- with(imp_state_year[[1]] %>% filter(ccode == 645, year == 2003),
                     c(latitude, longitude))
cf_iraq_bilat <- coalition_cf(dispute_dyads = dat_iraq_bilat,
                              fit_list = fit_base,
                              latlong = latlong_iraq,
                              comment = "USA only")

## Coalition of the willing
dat_iraq_base <- make_dispute_dyads(a = willing_coalition,
                                    b = c(645),
                                    year = 2003)
cf_iraq_base <- coalition_cf(dispute_dyads = dat_iraq_base,
                             fit_list = fit_base,
                             latlong = latlong_iraq,
                             comment = "+ Coalition of the willing")

## Baseline + key NATO allies
dat_iraq_nato <- make_dispute_dyads(a = c(willing_coalition, 220, 255, 640),
                                    b = c(645),
                                    year = 2003)
cf_iraq_nato <- coalition_cf(dispute_dyads = dat_iraq_nato,
                             fit_list = fit_base,
                             latlong = latlong_iraq,
                             comment = "+ Key NATO allies")

## All that + China and Russia
dat_iraq_total <- make_dispute_dyads(a = c(willing_coalition, 220, 255, 640, 365, 710),
                                     b = c(645),
                                     year = 2003)
cf_iraq_total <- coalition_cf(dispute_dyads = dat_iraq_total,
                              fit_list = fit_base,
                              latlong = latlong_iraq,
                              comment = "+ Russia and China")

dat_iraq <- collapse_results(cf_iraq_bilat, cf_iraq_base, cf_iraq_nato, cf_iraq_total)

## Plot force multipliers
pdat_iraq_fm <- dat_iraq$participant %>%
  left_join(country_names, by = "ccode") %>%
  group_by(country) %>%
  summarise(ccode = ccode[1], ratio = ratio[1]) %>%
  mutate(multiplier = 1 / ratio,
         category = case_when(
           ccode %in% c(2, 645) ~ 1,
           ccode %in% !!willing_coalition ~ 2,
           ccode %in% c(220, 255, 640) ~ 3,
           TRUE ~ 4
         ),
         category = factor(category,
                           levels = 1:4,
                           labels = c("Core players",
                                      "Coalition of the willing",
                                      "Key NATO allies",
                                      "UNSC vetoes")),
         country = recode(country, "Netherlands" = "Neth.")) %>%
  ungroup()

plot_iraq_fm <-
  ggplot(pdat_iraq_fm, aes(x = country, y = multiplier)) +
  geom_bar(stat = "identity") +
  facet_grid(~ category, scales = "free_x", space = "free") +
  scale_x_discrete("") +
  scale_y_continuous("Effective force multiplier")

pdf(file = "figures/figure_3.pdf",
    width = 7.5,
    height = 3.25)
print(plot_iraq_fm)
dev.off()

## Table of victory probabilities and bargaining quantities
tdat_iraq <- dat_iraq$dispute %>%
  mutate_at(vars(starts_with("pr")), ~ sprintf("%.2f", .)) %>%
  mutate_at(vars(starts_with("loc")), ~ sprintf("%.1f", .)) %>%
  mutate_at(vars(starts_with("scl")), ~ sprintf("%.1f", .)) %>%
  mutate_if(is.character, ~ str_replace(., "-", "$-$")) %>%
  mutate(entry = str_glue("{comment} & {pr_win_a} & {loc_a} & {scl_a} & {loc_b} & {scl_b} & {pr_war} \\\\")) %>%
  pull(entry) %>%
  c("\\begin{tabular}{lrrrrr|r}",
    "\\toprule",
    "Coalition & $p_A$ & $\\mu_A$ & $\\sigma_A$ & $\\mu_B$ & $\\sigma_B$ & $\\Pr$(War) \\\\",
    "\\midrule",
    .,
    "\\bottomrule",
    "\\end{tabular}")

writeLines(tdat_iraq, con = "tables/table_6.tex")


## -----------------------------------------------------------------------------
## World War II, European front (258)
## -----------------------------------------------------------------------------

dat_ww2_full <- make_dispute_dyads(a = c(2, 200, 220, 290, 365),
                                   b = c(255, 325),
                                   year = 1939)
cf_ww2_full <- coalition_cf(dispute_dyads = dat_ww2_full,
                            fit_list = fit_base,
                            mid_id = 258,
                            comment = "With both")

dat_ww2_init <- dat_ww2_full %>%
  filter(!(ccode1 %in% c(2, 365)),
         !(ccode2 %in% c(2, 365)))
cf_ww2_init <- coalition_cf(dispute_dyads = dat_ww2_init,
                            fit_list = fit_base,
                            mid_id = 258,
                            comment = "Baseline")

dat_ww2_usa <- dat_ww2_full %>%
  filter(ccode1 != 365,
         ccode2 != 365)
cf_ww2_usa <- coalition_cf(dispute_dyads = dat_ww2_usa,
                           fit_list = fit_base,
                           mid_id = 258,
                           comment = "With USA")

dat_ww2_rus <- dat_ww2_full %>%
  filter(ccode1 != 2,
         ccode2 != 2)
cf_ww2_rus <- coalition_cf(dispute_dyads = dat_ww2_rus,
                           fit_list = fit_base,
                           mid_id = 258,
                           comment = "With USSR")

dat_ww2 <- collapse_results(cf_ww2_init, cf_ww2_rus, cf_ww2_usa, cf_ww2_full)
dat_ww2_part <- dat_ww2[["participant"]]
dat_ww2_disp <- dat_ww2[["dispute"]]

pdat_ww2_part <-
  complete(dat_ww2_part, ccode, comment, fill = list(effort = 0, wt_effort = 0)) %>%
  left_join(country_names, by = "ccode") %>%
  mutate(alliance = if_else(ccode %in% c(255, 325), "Axis", "Allies"),
         Players = factor(comment, levels = c("Baseline", "With USSR", "With USA", "With both")),
         country = if_else(country == "Russia", "USSR", country))

plot_ww2 <- ggplot(pdat_ww2_part, aes(x = country, y = wt_effort)) +
  geom_col(aes(fill = Players),
           position = position_dodge2(preserve = "single")) +
  facet_grid(~ alliance, scales = "free_x", space = "free") +
  scale_x_discrete("") +
  scale_y_continuous("Equilibrium contribution") +
  scale_fill_manual(values = brewer.pal(5, "Reds")[-1]) +
  ggtitle("Counterfactuals for World War II, European theater, 1939")

pdf(file = "figures/figure_A2.pdf",
    width = 7,
    height = 3.25)
print(plot_ww2)
dev.off()

tab_ww2 <- make_table(dat_ww2_disp, win_name = "Allies Win")
writeLines(tab_ww2, con = "tables/table_A5.tex")


## -----------------------------------------------------------------------------
## World War II, Pacific front
## -----------------------------------------------------------------------------

## Specifying distance manually (centering in Japan) since no directly
## corresponding MID (WW2 MID gives the European front)
dat_pac_full <- make_dispute_dyads(a = c(2, 200, 365, 710),
                                   b = 740,
                                   year = 1941)
latlong_pac <- with(imp_state_year[[1]] %>% filter(ccode == 740, year == 1941),
                    c(latitude, longitude))
cf_pac_full <- coalition_cf(dispute_dyads = dat_pac_full,
                            fit_list = fit_base,
                            latlong = latlong_pac,
                            comment = "With USSR")

dat_pac_init <- dat_pac_full %>%
  filter(ccode1 != 365, ccode2 != 365)
cf_pac_init <- coalition_cf(dispute_dyads = dat_pac_init,
                            fit_list = fit_base,
                            latlong = latlong_pac,
                            comment = "Baseline")

dat_pac <- collapse_results(cf_pac_init, cf_pac_full)
dat_pac_part <- dat_pac[["participant"]]
dat_pac_disp <- dat_pac[["dispute"]]

tab_pac <- make_table(dat_pac_disp, win_name = "Allies Win")
writeLines(tab_pac, con = "tables/table_A6.tex")


## -----------------------------------------------------------------------------
## Taiwan Strait crisis (bargaining counterfactuals)
##
## MID #2987 Coding Taiwan as Side A, contrary to the MID coding, so that we can estimate
## the optimal offer by Taiwan
## -----------------------------------------------------------------------------

## Latitude and longitude for Taiwan Strait (via Wikipedia)
latlong_tws <- c(24.811111, 119.928333)

## Just China and Taiwan (which is how it's coded in MIDs)
dat_tws_bilat <- make_dispute_dyads(b = 710, a = 713, year = 1954)
cf_tws_bilat <- coalition_cf(dispute_dyads = dat_tws_bilat,
                             fit_list = fit_base,
                             latlong = latlong_tws,
                             comment = "Baseline")

## Only add US intervening in favor of Taiwan
dat_tws_us <- make_dispute_dyads(b = 710, a = c(2, 713), year = 1954)
cf_tws_us <- coalition_cf(dispute_dyads = dat_tws_us,
                          fit_list = fit_base,
                          latlong = latlong_tws,
                          comment = "With US")

## Only add USSR intervening in favor of China
dat_tws_ussr <- make_dispute_dyads(b = c(365, 710), a = 713, year = 1954)
cf_tws_ussr <- coalition_cf(dispute_dyads = dat_tws_ussr,
                             fit_list = fit_base,
                             latlong = latlong_tws,
                             comment = "With USSR")

## Add both major Cold War powers intervening
dat_tws_both <- make_dispute_dyads(b = c(365, 710), a = c(2, 713), year = 1954)
cf_tws_both <- coalition_cf(dispute_dyads = dat_tws_both,
                             fit_list = fit_base,
                             latlong = latlong_tws,
                             comment = "With both")

dat_tws <- collapse_results(cf_tws_bilat,
                            cf_tws_us,
                            cf_tws_ussr,
                            cf_tws_both)

## Plot force multipliers
pdat_tws_fm <- dat_tws$participant %>%
  left_join(country_names, by = "ccode") %>%
  mutate(country = if_else(country == "Russia", "USSR", country)) %>%
  group_by(country) %>%
  summarise(ccode = ccode[1], ratio = ratio[1], .groups = "drop") %>%
  mutate(multiplier = 1 / ratio,
         country = factor(country, levels = c("Taiwan", "China", "USA", "USSR")),
         category = fct_collapse(country,
                                 "Core Participants" = c("Taiwan", "China"),
                                 "Outside Allies" = c("USA", "USSR")))

plot_tws_fm <-
  ggplot(pdat_tws_fm, aes(x = country, y = multiplier)) +
  geom_bar(stat = "identity") +
  facet_grid(~ category, scales = "free_x") +
  scale_x_discrete("") +
  scale_y_continuous("Effective force multiplier")

pdf(file = "figures/figure_A3a.pdf",
    width = 4.5,
    height = 3.25)
print(plot_tws_fm)
dev.off()

## Table of victory probabilities and bargaining quantities
tdat_tws <- dat_tws$dispute %>%
  mutate_at(vars(starts_with("pr")), ~ sprintf("%.2f", .)) %>%
  mutate_at(vars(starts_with("loc")), ~ sprintf("%.1f", .)) %>%
  mutate_at(vars(starts_with("scl")), ~ sprintf("%.1f", .)) %>%
  mutate_at(vars(starts_with("x")), ~ sprintf("%.1f", .)) %>%
  mutate_if(is.character, ~ str_replace(., "-", "$-$")) %>%
  mutate(entry = str_glue("{comment} & {pr_win_a} & {x_star} & {loc_a} & {scl_a} & {loc_b} & {scl_b} & {pr_war} \\\\")) %>%
  pull(entry) %>%
  c("\\begin{tabular}{lrrrrrr|r}",
    "\\toprule",
    "Participants & $p_A$ & $x^*$ & $\\mu_A$ & $\\sigma_A$ & $\\mu_B$ & $\\sigma_B$ & $\\Pr$(War) \\\\",
    "\\midrule",
    .,
    "\\bottomrule",
    "\\end{tabular}")

writeLines(tdat_tws, con = "tables/table_A7a.tex")

## Hypothetical "present day" (2012 is last year in our data) version
dat_twsnow_bilat <- make_dispute_dyads(b = 710, a = 713, year = 2012)
cf_twsnow_bilat <- coalition_cf(dispute_dyads = dat_twsnow_bilat,
                             fit_list = fit_base,
                             latlong = latlong_tws,
                             comment = "Baseline")
dat_twsnow_us <- make_dispute_dyads(b = 710, a = c(2, 713), year = 2012)
cf_twsnow_us <- coalition_cf(dispute_dyads = dat_twsnow_us,
                          fit_list = fit_base,
                          latlong = latlong_tws,
                          comment = "With US")
dat_twsnow_ussr <- make_dispute_dyads(b = c(365, 710), a = 713, year = 2012)
cf_twsnow_ussr <- coalition_cf(dispute_dyads = dat_twsnow_ussr,
                             fit_list = fit_base,
                             latlong = latlong_tws,
                             comment = "With Russia")
dat_twsnow_both <- make_dispute_dyads(b = c(365, 710), a = c(2, 713), year = 2012)
cf_twsnow_both <- coalition_cf(dispute_dyads = dat_twsnow_both,
                             fit_list = fit_base,
                             latlong = latlong_tws,
                             comment = "With both")
dat_twsnow <- collapse_results(cf_twsnow_bilat,
                               cf_twsnow_us,
                               cf_twsnow_ussr,
                               cf_twsnow_both)

## Plot force multipliers
pdat_twsnow_fm <- dat_twsnow$participant %>%
  left_join(country_names, by = "ccode") %>%
  group_by(country) %>%
  summarise(ccode = ccode[1], ratio = ratio[1], .groups = "drop") %>%
  mutate(multiplier = 1 / ratio,
         country = factor(country, levels = c("Taiwan", "China", "USA", "Russia")),
         category = fct_collapse(country,
                                 "Core Participants" = c("Taiwan", "China"),
                                 "Outside Allies" = c("USA", "Russia")))

plot_twsnow_fm <-
  ggplot(pdat_twsnow_fm, aes(x = country, y = multiplier)) +
  geom_bar(stat = "identity") +
  facet_grid(~ category, scales = "free_x") +
  scale_x_discrete("") +
  scale_y_continuous("Effective force multiplier")

pdf(file = "figures/figure_A3b.pdf",
    width = 4.5,
    height = 3.25)
print(plot_twsnow_fm)
dev.off()

## Table of victory probabilities and bargaining quantities
tdat_twsnow <- dat_twsnow$dispute %>%
  mutate_at(vars(starts_with("pr")), ~ sprintf("%.2f", .)) %>%
  mutate_at(vars(starts_with("loc")), ~ sprintf("%.1f", .)) %>%
  mutate_at(vars(starts_with("scl")), ~ sprintf("%.1f", .)) %>%
  mutate_at(vars(starts_with("x")), ~ sprintf("%.1f", .)) %>%
  mutate_if(is.character, ~ str_replace(., "-", "$-$")) %>%
  mutate(entry = str_glue("{comment} & {pr_win_a} & {x_star} & {loc_a} & {scl_a} & {loc_b} & {scl_b} & {pr_war} \\\\")) %>%
  pull(entry) %>%
  c("\\begin{tabular}{lrrrrrr|r}",
    "\\toprule",
    "Participants & $p_A$ & $x^*$ & $\\mu_A$ & $\\sigma_A$ & $\\mu_B$ & $\\sigma_B$ & $\\Pr$(War) \\\\",
    "\\midrule",
    .,
    "\\bottomrule",
    "\\end{tabular}")

writeLines(tdat_twsnow, con = "tables/table_A7b.tex")


date()
